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Abstract 

Transport properties of liquid methanol and ethanol are predicted by mole- 
cular dynamics simulation. The molecular models for the alcohols are rigid, 
non-polarizable and of united-atom type. They were developed in preced- 
ing work using experimental vapor-liquid equilibrium data only. Self- and 
Maxwell-Stefan diffusion coefficients as well as the shear viscosity of metha- 
nol, ethanol and their binary mixture are determined using equilibrium molec- 
ular dynamics and the Green-Kubo formalism. Non-equilibrium molecular 
dynamics is used for predicting the thermal conductivity of the two pure 
substances. The transport properties of the fluids are calculated over a wide 
temperature range at ambient pressure and compared with experimental and 
simulation data from the literature. Overall, a very good agreement with the 
experiment is found. For instance, the self-diffusion coefficient and the shear 
viscosity are predicted with average deviations of less 8% for the pure alco- 
hols and 12% for the mixture. The predicted thermal conductivity agrees 
in average within 5% with the experimental data. Additionally, some ve- 
locity and shear viscosity autocorrelation functions are presented and dis- 
cussed. Radial distribution functions for ethanol are also presented. The 
predicted excess volume, excess enthalpy and the vapor-liquid equiUbrium 
of the binary mixture methanol + ethanol are assessed and the vapor-liquid 
equilibrium agree well with experimental data. 

Keywords: Green-Kubo; reverse NEMD; self-diffusion coefficient; Maxwell- 
Stefan diffusion coefficient; shear viscosity; thermal conductivity. 
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1 Introduction 



There has been significant progress in recent years to study and predict both static 
and dynamic thermophysical properties of fluids by molecular simulation. This is 
particularly interesting for transport properties of liquids, where, due to the com- 
plexity of the involved physical mechanisms, other theoretical approaches are of- 
ten unsatisfactory. Hence, molecular dynamics is a promising method for predict- 
ing transport properties of liquids. However, the number of publications regarding 
transport properties for technically relevant liquid systems that achieve quantita- 
tively accurate results is still low. The majority of these publications deals with the 
self-diffusion coefficient, leaving aside the more demanding Maxwell-Stefan dif- 
fusion coefficient, shear viscosity and thermal conductivity. Furthermore, molec- 
ular simulation studies on transport properties for complex liquids, e.g. mixtures 
containing highly polar hydrogen-bonding molecules, are rare. This is also the 
case for methanol, ethanol and their binary mixture. It should be pointed out 
that there is a series of articles on transport properties of pure liquid ethanol by 
Petravic et al. [|Tl|2l[3l|. There are several works on transport properties of associ- 
ating mixtures: Typically, aqueous solutions of methanol [IH [51 [6l |7l [H or ethanol 
[|7l [H m [lOl [m have been studied using molecular dynamics simulation. How- 
ever, to our knowledge the mixture methanol -i- ethanol has not yet been regarded 
in that sense. 

Two fundamentally different methods for calculating transport properties by 
molecular dynamics simulation are available. The equilibrium methods (EMD), 
using either the Green-Kubo formalism or the Einstein relations, analyze the time 
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dependent response of a fluid system to spontaneous fluctuations. With non- 
equilibrium molecular dynamics (NEMD), on the other hand, the system's re- 
sponse to an externally applied perturbation is analyzed. NEMD was developed 
in order to increase the signal to noise ratio and to improve statistics and con- 
vergence. Both methods have different advantages and disadvantages, but are 
comparable in efficiency, as shown e.g. by Dysthe et al. [fTlll . 

The success of molecular dynamics simulation to predict thermodynamic prop- 
erties is highly determined by the potential model used to describe the intermolec- 
ular interactions. Numerous molecular models for methanol and ethanol have 
been proposed in the literature [l[Il[Il[l5l[ia[I71[Il[iaiMl2DIMl23llMl25l. 
However, the ability of these models for predicting transport properties has been 
probed predominantly regarding the self-diffusion coefficient, cf. Table [T] 

For methanol, mainly the rigid OPLS molecular model from Jorgensen [fTSlI 
and its modifications proposed by Haughney et al. [fT4l as well as by van Leeuwen 
and Smit fJ6J . have been used to predict the self-diffusion coefficient. It was 
shown by comparison of the models by Jorgensen [fT3l and by Haughney et al. 
[[T4]| that the rigid OPLS model and its modifications perform well for the self- 
diffusion coefficient, deviations to experimental data range from 7 to 20% [fT4l . 
Van de Ven-Lucassen et al. O showed that the model from van Leeuwen and 
Smit [fT6l predicts the self-diffusion coefficient of methanol at ambient tempera- 
ture more accurately than the polarizable model of Caldwell and KoUman [fTSl . 
and than the models by Haughney et al. [fT4ll . Wheeler et al. [26] predicted 
the shear viscosity of methanol within 10% of experimental values. The thermal 
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conductivity of methanol at ambient conditions was reported by Nieto-Draghi Q 
with a deviation of + 1 1 %. 

Different rigid |[IS|271, semi-flexible [28], flexible ^ and polarizable CH 
l29l molecular models for ethanol have been used to predict the self-diffusion co- 
efficient. Thereby, it has been found that rigid, united-atom molecular models pre- 
dict the self-diffusion coefficient more reliably than polarizable, flexible models. 
Most ethanol models assessed by Wang and Cann ll25l as well as the Jorgensen 
model [fT3l investigated by Saiz et al. Il28l overestimate the self-diffusion coef- 
ficient rather significantly with deviations ranging from -i-8 to -1-50%. Allowing 
a potential model for bending results in a further increase in the predicted val- 
ues Gonzalez et al. Il30l found that rigid, non-polarizable models perform 
better than the polarizable ones for predictions of the self-diffusion coefficient. 
The shear viscosity of ethanol has been predicted using EMD [|2l and NEMD 
[|3T1l . Both works achieved good results (average deviations are 11 and 7%, re- 
spectively) using the united-atom model and the all-atom model of Jorgensen et 
al. [fT3l [T9l , respectively. Simulations performed by Petravic [3 J underestimate 
the thermal conductivity of liquid ethanol near ambient temperatures by approxi- 
mately -10%. 

The present article assesses the predictive power of rigid, non-polarizable, 
united-atom type Lenard- Jones based models with superimposed point charges 
for methanol and ethanol developed in prior work of our group [|32l [33l . These 
models have been optimized on the basis of experimental data of vapor pressure 
and saturated liquid density in the whole temperature range from the triple point 
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to the critical point. No experimental data for transport properties was thereby 
taken into account. The most significant transport properties, i.e. self-diffusion 
coefficient, shear viscosity and thermal conductivity, were regarded here for the 
pure liquids using EMD and NEMD. Furthermore, the two molecular models were 
tested regarding the prediction of self- and Maxwell-Stefan diffusion coefficients 
and shear viscosity of the binary mixture. Additionally, binary vapor-liquid equi- 
libria as well as excess volume and excess enthalpy were predicted. The goal of 
this study is to demonstrate the ability of rigid, non-polarizable molecular models, 
adjusted to pure substance vapor-liquid equilibria only,to accurately describe the 
transport properties of strongly polar and H-bonding liquids. 

This work is organized as follows: Firstly, the employed molecular models 
are briefly described. Then, the simulation techniques are explained, the results 
for self -diffusion coefficient, shear viscosity and thermal conductivity of the pure 
liquids are presented and compared to experimental and simulation data from the 
literature. The self- and Maxwell-Stefan diffusion coefficients as well as the shear 
viscosity for the mixture methanol + ethanol are presented and compared with 
experimental data where possible. Furthermore, the velocity and shear viscosity 
autocorrelation functions are described and discussed. Simulated radial distribu- 
tion functions of ethanol are also presented. Thereafter, the simulation results for 
vapor-liquid equilibria, excess volume and excess enthalpy for the binary mixture 
are given in comparison to experimental values. Finally, conclusions are drawn. 
Simulation details are summarized in the Appendix. 
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2 Molecular models 



Throughout this work, rigid, non-polarizable molecular models of united-atom 
type from earlier work of our group [|32l[33|| were used. Both models account for 
the intermolecular interactions, including hydrogen bonding, by a set of united- 
atom Lennard- Jones (LJ) sites and superimposed point charges. Hence, the po- 
tential energy Uij between two alcohol molecules i and j can be written as 



where a is the site index of molecule b the site index of molecule j, and n 
and m are the numbers of interaction sites of molecules i and j, respectively, rijab 
represents the site-site distance between molecules / and j. The LJ size and energy 
parameters are Oab and Eab- lia and qjt are the point charges located at the sites a 
and b on the molecules / and j, and Eq is the permittivity of vacuum. 

The methanol model has two LJ sites, i.e. one for the methyl group and one 
for the hydroxyl group. Additionally, three point charges are superimposed, i.e. 
one at each of the LJ site centers and one at the nucleus position of the hydroxyl 
hydrogen. Ethanol was modeled by three LJ sites, i.e. one for the methyl, the 
methylene and the hydroxyl group each. Three point charges are located here as 
well on the methylene and hydroxyl LJ site centers and on the nucleus position of 
the hydroxyl hydrogen. 

The parameter set of the methanol model was obtained using the optimiza- 
tion procedure of StoU ll34l on the basis of the geometry of the model from van 




(1) 
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Leeuwen and Smit lfT6l . The geometry of the ethanol model was derived from 
ab initio quantum chemistry calculations. The AUA4 parameters of Ungerer et 
al. [|35l were applied for the methyl and methylene LJ sites. The point charge, 
LJ size and energy parameters as well as the offset of the hydroxyl LJ site were 
fitted to experimental vapor pressure and saturated liquid density. It should be 
pointed out that no information on transport or mixture properties was included 
in the optimization procedures of any molecular model employed here so that the 
data presented below is strictly predictive. The model parameters were taken from 
[l32l[33ll and are summarized in Table [21 

To define a molecular model for a binary mixture on the basis of pairwise ad- 
ditive pure substance models, only the unlike interactions have to be specified. In 
case of polar interaction sites, i.e. point charges here, this can straightforwardly be 
done using the laws of electrostatics. However, for the unlike LJ parameters there 
is no physically sound approach [|36ll . and for predictions combination rules have 
to be employed. Here, the interaction between unlike LJ sites of two molecules 
are calculated applying the Lorentz-Berthelot combining rules: 
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and 
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3 Methodology 



Dynamic properties can be obtained from EMD simulations by means of the 
Green-Kubo formalism [|37l [38l . These equations are a direct relationship be- 
tween a transport coefficient and the time integral of an autocorrelation function 
of a particular microscopic flux in a system in equilibrium. This method was used 
in this work to calculate self- and Maxwell-Stefan diffusion coefficients as well as 
the shear viscosity. The thermal conductivity was calculated using a modified ver- 
sion of the boundary driven (BD) reverse NEMD algorithm from Miiller-Plathe 
[|39ll . also referred to as PeX algorithm. 

3.1 Diffusion coefficients 

The self-diffusion coefficient is related to the mass flux of single molecules 
within a fluid. Therefore, the relevant Green-Kubo expression is based on the 
individual molecule velocity autocorrelation function as follows 



where Vyt(?) is the center of mass velocity vector of molecule k at some time t, and 
<...> denotes the ensemble average. Eq. dU is an average over all molecules 
in a simulation, since all contribute to the self-diffusion coefficient. In a binary 
mixture, the Maxwell- Stefan diffusion coefficient D/y can also be written in terms 
of the center of mass velocity 




(4) 
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where Mi and Xi are the molar mass and mole fraction of species i. From the 
expression above, the collective character of the Maxwell- Stefan diffusion coeffi- 
cient is evident. 

Since the present simulations provide self- as well as the Maxwell-Stefan dif- 
fusion coefficients simultaneously, a comparison to the simple predictive approach 
suggested by Darken Il40l is possible. The Maxwell-Stefan diffusion coefficient 
Dij is estimated from the self-diffusion coefficients of the two components and 
Dj in the binary mixture ll40l 

Dij=XiDi + XjDj- (6) 

This approach was found to be superior to the correlations by Vignes [[4T]| and 
by Caldwell and Babb [|42| in prior work on three binary mixtures B3l . How- 
ever, it is of little practical use due to the fact that it relies on the self-diffusion 
coefficients in the mixture, which are rarely available. 

3.2 Shear viscosity 

The shear viscosity r\, as defined by Newton's "law" of viscosity, is a measure of 
the resistance of a fluid to a shearing force p4|. It is associated with momentum 
transport under the influence of velocity gradients. Hence, the shear viscosity can 
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be related to the time autocorrelation function of the off-diagonal elements of the 
stress tensor [|45l 

^ = ^/^dr (7,^(0 -/^(O)), (7) 

where V stands for the molar volume, kg is the Boltzmann constant, and T denotes 
the temperature. Averaging over all three independent elements of the stress ten- 
sor, i.e. 7^, 7p^,and J^p, improves the statistics of the simulation. The component 
of the microscopic stress tensor Jp is given in terms of the molecule positions 
and velocities by ll44l 



i=l ^ i=\ j^ik=\l=\ "'kl 

Here, the lower indices / and k count the interaction sites, and the upper indices 
X and y denote the spatial vector components, e.g. for velocity vf or site-site dis- 
tance r^j. The first term of Eq. ([8]) is the kinetic energy contribution and the 
second term is the potential energy contribution to the shear viscosity. Conse- 
quently, the Green-Kubo integral (Q can be decomposed into three parts, i.e. 
the kinetic energy contribution, the potential energy contribution and the mixed 
kinetic-potential energy contribution [|44|. Eqs. ^ and ([8]) may directly be ap- 
plied to mixtures. 
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3.3 Thermal conductivity 

The thermal conductivity X, as defined by Fourier's "law" of heat conduction, 
characterizes the capability of a substance for molecular energy transport driven 
by temperature gradients. In BD-NEMD simulation used in this work, two slabs of 
a given thickness, large enough to contain a sufficiently large number of molecules 
on average but much smaller than the total length of the simulation volume in z 
direction, are defined. In order to create a heat flux in z direction, the molecule 
with the highest kinetic energy in the cold slab is selected with a frequency D 
to perform a hypothetical elastic collision with the molecule having the lowest 
kinetic energy in the hot slab ll46l . After such a virtual collision the two molecules 
have exchanged their momentum without changing their respective position in 
space. The new velocity of the molecule in the cold region is then 



and, for the molecule in the hot region 



^.^^_,..^ '>''<" + '>".<\ (10) 

mc + nih 

where nic and w/, are the respective masses of the selected molecules in the cold 
and hot slabs, v'^^'^ and v^^"' stand, respectively, for the velocities of the molecules 
before and after the collision. 

During the elastic collision the molecules exchange their momentum. Thus, 
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the total momentum and the total energy of the system remain constant. The ki- 
netic energy exchanged by means of these elastic collisions determines the energy 
flux in the system. The heat flow density in the steady state is then given by 



^''^^^ = 'lAt ^ "^[^^rf-ivff)- (11) 

transfers 

In this expression, A is the cross sectional area in x and y direction, t is the time 
interval of measure in the steady state during which a given number of momentum 
transfers have occurred. {Jc)f is the heat flux in z direction, which can indistinctly 
be evaluated in the cold or the hot slab. It is required that the number of molecules 
in both slabs is sufficiently large so that the intermolecular interactions within the 
slabs will properly thermalize the velocity distribution before a new collision takes 
place. The thermal gradient Vr^ is obtained from a linear fit of the temperature 
profiles from the simulation. In the steady state, the thermal conductivity can be 
obtained by 



4 Simulation results 
4.1 Diffusion coefficients 

The self-diffusion coefficient of pure methanol and pure ethanol was predicted at 
ambient pressure in the temperature range from 210 to 340 K. Present numerical 
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data are given in Table [3l Figures \T\ and [2] show the self-diffusion coefficient for 
methanol and ethanol as a function of temperature, respectively. The statistical 
error of the simulation data is estimated to be in the order of 1%. 

The self-diffusion coefficient of methanol shows a very good agreement with 
the experimental values. At temperatures below 280 K, the self-diffusion coef- 
ficient was underpredicted by approximately —7%. Above ambient temperature, 
the mobility of the model molecules increases slightly more rapidly than that of 
real methanol molecules. Hence, for temperatures above 300 K the self-diffusion 
coefficient was overestimated on average by +8%. However, the significant scat- 
ter of the experimental data should be taken into account. In comparison to pre- 
vious simulation studies from the literature, the present results cover the widest 
temperature range. The broadest available study by Haughney et al. ffT4l . using a 
modified Jorgensen model, is comparable with present results near ambient tem- 
perature where average deviations are 9 %. However, the data from Haughney et 
al. is scattering and the self-diffusion coefficient is significantly underestimated 
above 300 K. 

There is an excellent agreement between experimental data and present self- 
diffusion coefficient for ethanol. The predicted data also correctly reproduces 
the temperature dependence over the whole temperature range from 239 to 338 
K. Contrary to previous simulation reports [[Il|25l|27l|28ll29l[30l|471, where the 
self-diffusion coefficient of ethanol has mostly been significantly overestimated, 
deviations of the present data are negative and on average only —6%. This is 
explained by the different potential model used in this work. 
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Self-diffusion coefficients of methanol and ethanol in their binary mixture 
were predicted at ambient conditions for the entire composition range. Numerical 
results are presented in Tableland plotted in Figure [3l The estimated statistical 
uncertainty of the simulation data is also 1%. As shown in Figure [3l there is an 
excellent agreement between predicted self-diffusion coefficients and the experi- 
mental data for both alcohols. Present simulation results slightly overestimate the 
self-diffusion of both alcohols by +5% on average. Furthermore, the predicted 
data reproduce the composition dependence well. 

The Maxwell-Stefan diffusion coefficient of the mixture methanol -i- ethanol 
calculated by molecular simulation is reported in Table |4l In Figure |4l the pre- 
dicted values are plotted as a function of the composition. The statistical uncer- 
tainties are on average 15%. These large error bars are attributed to the collective 
nature of the Maxwell-Stefan diffusion coefficient. Due to the lack of experi- 
mental data, present simulation results are only compared to the predictions of 
Darken's model HOll . cf. Eq. I©. This model can be considered to be a good 
approximation here, in agreement to previous findings for three other mixtures 

m. 

4.2 Shear viscosity 

The shear viscosity of pure methanol and pure ethanol was predicted at ambient 
pressure for a temperature range of about 100 K around ambient conditions and is 
reported in Table [31 The temperature dependence of the shear viscosity is shown 
in Figures [5] and [6] for methanol and ethanol, respectively. Statistical errors are in 
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the range from 10 to 17%. 

In the case of methanol, the predicted shear viscosity agrees very well with 
the experimental data over the whole regarded temperature range with an average 
deviation of 8%. At higher temperatures, the predicted viscosities are closer to the 
experimental ones than at lower temperatures, cf. Figure [5l It should be pointed 
out that the highest deviation, being 20%, was found at the lowest regarded tem- 
perature 220 K. Here, the high density in combination with strongly interacting 
molecules causes large simulation uncertainties due to the long-time behavior of 
the shear viscosity autocorrelation function. Present shear viscosity predictions 
are better than the EMD data of Wheeler et al. Il26l and are comparable with the 
NEMD Edwald data of the same author. 

The calculated shear viscosity for ethanol shows a good agreement to exper- 
imental data with an average deviation of 8%. It can be noticed in Figure |6] that 
there is a tendency to underestimate the shear viscosity. Petravic and Delhom- 
melle [El also underestimated the shear viscosity using the OPLS model from 
Jorgensen fT3\. At high densities, however, present data better predict the shear 
viscosity. Zhao et al. OTI recently reported shear viscosity data using the OPLS- 
AA model [fT9l and NEMD which is in very good agreement with experimental 
data. 

The three contributions to the shear viscosity were also calculated. In agree- 
ment with the results of Meier et al. [|48ll . the potential energy contribution is 
dominant for both studied alcohols, r[pp > 90%. In general, for the studied tem- 
perature range, the kinetic energy term has a negligible positive contribution to 
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the shear viscosity, r[kk < 2%. The mixed kinetic-potential energy term may con- 
tribute positively or negatively to the shear viscosity. Here, this contribution varies 
between 0.1 and 8%. 

The simulated shear viscosity data of the mixture methanol -i- ethanol is given 
in Table |4] and plotted in Figure U\ together with the available experimental data. 
There is a very good agreement between experimental and predicted shear vis- 
cosity, being mostly within the simulation uncertainty of around 12%. Deviations 
from the experimental data are on average 12% as well. 

4.3 Thermal conductivity 

Present results for the thermal conductivity of pure methanol and pure ethanol at 
ambient pressure are summarized in Tabled The agreement with the experimen- 
tal data is remarkable with an average deviation of 5% for methanol and of 2% 
for ethanol, cf. Figures [8] and [9l The thermal conductivity for methanol is slightly 
underestimated by the present predictions, while for ethanol it is slightly overesti- 
mated. Previous results for the thermal conductivity of methanol Q and ethanol 
[|3l|, using the OPLS models, show deviations of about 11 and 10%, respectively. 
In contrast to present results, the calculated thermal conductivity by Petravic [[3]| 
lies below the experimental data. 
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4.4 Time dependent behavior of the Green-Kubo integrands 
Velocity autocorrelation function 

The velocity autocorrelation function provides an interesting insight into the liq- 
uid state. Selected normalized velocity autocorrelation functions (VACF) of pure 
methanol and ethanol are shown in Figure \Wi In agreement with findings in the 
literature ^49\. VACF may have two minima at short times (f < 0.5 ps). At low 
temperatures and high densities, the VACF of both pure liquids decrease sharply 
and become negative. After this first minimum, they increase slightly to drop 
again into a lower minimum, for methanol, or to a higher minimum, for ethanol. 
Beyond the second minimum, the VACF converge to zero following a character- 
istic path. The observed negative values of the VACF are related to backscat- 
tering collisions, also known as cage-effect, which govern short time dynamics. 
Here, the free movement of the molecules is hindered by surrounding molecules, 
which form a cage-like structure. Michels and Trappeniers [|50l attributed this 
phenomenon to the formation of "bound states" since this effect can only be ob- 
served for molecular models that include attractive forces. 

For higher temperatures at ambient pressure, where the liquid density is lower, 
molecular collisions are less frequent and the depth of the minima is gradually 
reduced until they almost disappear, cf. Figure[lOl At short times, the VACF decay 
faster at lower temperatures. Moreover, the initial decay of the VACF of methanol 
occurs earlier than in the corresponding VACF of ethanol. Methanol VACF exhibit 
a more pronounced backscattering oscillations than those of ethanol, particularly 
at low temperatures near its melting temperature (r„, = 175.37 K) [1511 . 
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Figure nn shows the behavior of the integrals of the VACF given by Eq. dH). 
It can be seen that the self-diffusion coefficient converges in all regarded cases to 
its final value after less than 4 ps. Thus, after this time interval the long-time tail 
of the VACF does not contribute significantly anymore. 

The behavior of the VACF of methanol and ethanol was also considered in the 
mixture. From Figures [121 and [13] it can be seen that the VACF of methanol and 
ethanol in the mixture differ from those of the pure alcohols at the same tempera- 
ture and pressure. Methanol VACF exhibit a significantly increased backscattering 
due to the presence of ethanol, cf . Figure [I2l This indicates that at short times, 
ethanol molecules further restrain the mobility of methanol molecules. Ethanol 
VACF also show a composition dependent behavior. In this case, the backscat- 
tering is reduced as the concentration of methanol increases, cf . Figure \T3\ This 
suggests that methanol molecules prevent the formation of cage-like structures 
here, allowing ethanol molecules to move more freely in the mixture. 

Shear viscosity autocorrelation function 

Alder et al. |[52l were the first to suggest the existence of a long-time tail in 
the shear viscosity autocorrelation function (SACF) near the phase boundary to 
solidification due to the low compressibility of such liquid states. Luo and Ho- 
heisel ll53l [54l showed that the long-time behavior of the SACF is determined 
primarily by long-lived orientational correlations, occurring in liquids contain- 
ing anisotropic molecules. This behavior causes convergence problems in the 
Green-Kubo integral, since the SACF must be calculated over a relatively long 
time interval. 
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Figure [141 shows selected normalized SACF of pure methanol. The presence 
of a significant long-time tail can clearly be noticed in Figure \T5\ where the inte- 
grals of the SACF of methanol are plotted for three temperatures. As expected, 
the behavior of the long-time tail of the SACF is highly dependent on molecule 
species and state point. For pure methanol at a low temperature of 260 K, approx- 
imately 60% of the contribution to the total shear viscosity comes from the first 
2 ps. At a higher temperature of 340 K the first 2 ps contribute 90% to the total 
shear viscosity. In the case of ethanol (not shown here), the SACF are significant 
for even longer times, i.e. at 263 K the first 2 ps of the SACF contribute only 45% 
to the final value of the shear viscosity. At 338 K this contribution reaches 75%. 
This difference can be explained by the higher molecular anisotropy of ethanol. 

The kinetic, potential and mixed kinetic-potential energy contributions to the 
SACF were also considered in this work. Over the whole regarded range of state 
points, the behavior of the kinetic contribution is analogous to the monotonic 
chair-like decay of autocorrelation functions for non- attractive fluids. The cor- 
responding time integral exhibits no important contribution from the long-time 
behavior, cf. Figure The potential energy contribution shows at short times 
oscillations similar to those of the VACF, cf. Figures[T0]and[l71 These oscillations 
have been also related to the presence of "bound states" by Michels and Trappe- 
niers [|50l and by Meier et al. [48 J. The long-time tail is significant and approxi- 
mates the behavior of the total SACF. The mixed kinetic -potential energy contri- 
bution shows a completely different behavior at short times, cf. Figure [HI where 
data is presented for pure methanol at four temperatures. At very short times 
(t < 0.1 ps), two roughly mirrored patterns can be seen, the SACF decay /raise 
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sharply to reach a minimum/maximum and oscillate strongly during their decay 
to zero. Meier et al. [|48l . who investigated simple LJ fluids, only observed initial 
positive peaks in the mixed kinetic-potential energy contribution. To our knowl- 
edge, the presence of negative initial peaks has not been reported before. The 
long-time tail behavior of the mixed contribution is difficult to analyze because of 
the presence of strong oscillations. 

4.5 Hydrogen bonding dynamics 

Petravic et al. (T) have shown that internal degrees of freedom play a crucial 
role on the relaxation of the hydrogen bonding dynamics of liquid ethanol. For 
instance, they found that "freezing" bending and torsion in the OPLS model of 
ethanol produces a slowdown of the dynamics of the hydrogen bonds at 273 K 
(they observed an increase of the hydrogen bond life time of about 27%). For the 
methanol model used in this work, an excellent agreement between simulated and 
experimental site-site radial distribution functions and hydrogen bond statistics 
has been reported in previous work [l32ll . In the present work, the hydrogen bond 
dynamics of the rigid ethanol model was compared with the results using a flexi- 
ble model version including bending and torsion. For this porpouse, a geometrical 
hydrogen bonding criterion with roH below 2.6 A, roo below 3.5 A and 0(^oH...o) 
below 30° together with a life time probability accounting for intermittent hy- 
drogen bonds [[55] was used. In the flexible version, bending and torsion angles 
were allowed to evolve using a standard harmonic and a Fourier series potential, 
respectively. The potential constants were taken from the TraPPE intermolecu- 
lar potential for alcohols [|22|. keeping the equilibrium values of the bending and 
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torsion angles of the rigid model and all the remaining inter-molecular model pa- 
rameters as well as the bond distances between all sites. Simulations were then 
performed at a specified temperature of 273 K (using the Nose-Hoover thermo- 
stat) and a density of 17.62 mol/1 to compute the radial distribution functions, the 
hydrogen bond dynamics and self-diffusion coefficient for both, rigid and flexible 
version of the present ethanol model. 

The comparison of the radial distribution functions from the flexible and the 
rigid model can be seen in Figure [191 It was found that there are only minor differ- 
ences in the liquid structure as revealed by the three radial distribution functions 
goo, gOH and gHH- The position and height of the first peaks are almost identi- 
cal for both models, with the exception of gHH, where a more pronounced peak 
is observed for the rigid model, cf. Figure \T9\ Only small differences are noted 
around the second peak of goo and gon, niainly related to the molecules hydrogen 
bonded in the second solvation shell. In the latter case, a more pronounced second 
peak was found for the rigid model. In general, it can be said that the inclusion 
of the flexibility does not strongly affect the structure of liquid ethanol around the 
first solvation shell. 

In the next step, differences in the hydrogen bond dynamics for the two mod- 
els were analyzed. Average life times of hydrogen bonds of 137 and 80 ps were 
obtained for the rigid and flexible versions, respectively. The flexibility signifi- 
cantly reduces the hydrogen bond life time by a factor of 1 .7. Petravic et al. found 
in a similar investigation for the OPLS ethanol model a reduction in the hydrogen 
bond dynamics of about 1.3. Hence, the effect is more pronounced for the present 
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model. The hydrogen bond dynamics has a significant impact on the transport 
properties. Eg. at 273 K a self-diffusion coefficient of 0.802 x lO^^m^/s was 
obtained for the flexible version, which is in perfect agreement with the values re- 
ported by Petravic et al. [IJ. This value is higher by a factor of 1.6 compared to the 
self-diffusion coefficient obtained for the rigid model (0.493 x lO^^m^/s) which 
is in excellent agreement with the experimental data (0.48 x lO^^m^/ s) ll56l . This 
comparison shows the impact of the internal degrees of freedom on the hydrogen 
bond dynamics and the self-diffusion coefficient. Similar effects can be expected 
also for other transport properties. Flexible and rigid models of H-Bonding sys- 
tems obviously show significant differences regarding transport properties. This 
should also be kept in mind when assessing the comparisons of different models 
shown throughout this work. 

4.6 Vapor-liquid equilibria 

Vapor-liquid equilibria (VLE) often serve as a benchmark to assess molecular 
models. This is a useful choice, due to the fact that VLE represent the dominant 
feature in the region of fluid states. Indeed, they are nowadays preferably being 
used in the development of molecular models by optimizing the model parameters 
to VLE properties. This approach was also followed in the development of both 
regarded pure fluid models. In the next step, molecular models can be assessed 
on how they perform in mixtures. Based on experience for numerous binaries and 
ternaries, very good success can be expected, if the pure substance models are 
good. 
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Calculated VLE data of pure methanol and ethanol using present potential 
models have been reported previously for a temperature range from 270 to 490 K 
[|32l[33l . The average deviations found for methanol in vapor pressure, saturated 
liquid density and heat of vaporization are of 1.1, 0.6, and 5.5%, respectively 
[|32ll . The corresponding reported average deviations for ethanol are of 3.7, 0.3 
and 0.9%, respectively ll33l . 

Present simulation results for the binary VLE are given in Table |5] and com- 
pared in Figures|20]and|2T|to experimental data by Butcher and Robinson [|57l and 
to the Peng-Robinson equation of state [|58l . As in the molecular simulations no 
adjustment to binary data was carried out, also the Peng-Robinson equation was 
used without any adjustment to such data (kfj = 0). Figure [201 shows the almost 
"ideal" behavior of methanol -i- ethanol for two temperatures. In fact, these were 
the highest two isotherms for which experimental data were available to us. At the 
higher one (443.15 K), however, it can be seen from the experimental data that the 
bubble line is not fully linear. The comparison of simulation and experiment is ex- 
cellent throughout, both the pressure slope and the composition relation between 
vapor and liquid match good, only minor deviations are found in the ethanol-rich 
composition range where slightly too high vapor compositions are obtained. 

4.7 Excess properties 

Excess properties are basic mixing properties which contain real mixing behavior. 
In case of methanol -i- ethanol non-ideal mixing effects are very weak. At equimo- 
lar composition, where the mixing influence is strongest, experimental data sug- 
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gest a volume increase of 0.017% only. Also for enthalpy, a tiny positive excess 
contribution of 0.011 % was found by experiment. This resolution is difficult to 
achieve by molecular simulation. 

To determine excess properties, a straightforward approach was used here. 
Three simulations at specified temperature and pressure were performed, two for 
the pure substances and one of the mixture at a given composition. A binary 
excess property given by 

= y-xi-yi-X2-y2, (13) 

where y^ represents the excess volume or the excess enthalpy h^. Enthalpies 
or volumes of the pure liquids methanol and ethanol as well as of their binary 
mixture are denoted by ji, j2 and y, respectively. 

Table [6] contains present simulation data for ambient conditions which are 
compared to experimental results [|59l [60l [6T1 [62ll in Figures [22l and [231 The 
excess volume, cf . Figure [22l agrees with the experiment in all cases within its 
seemingly large error bars. However, it should be noted that the simulation uncer- 
tainty is only around 0.1% of the total volume of the mixture. A slight tendency of 
the simulation data towards higher excess volumes might be assumed. In Figure 
[23l the present excess enthalpy is shown together with experimental data. The 
maximal excess contribution from experiment is only 0.01 1%, which is well be- 
low the simulation uncertainty of around 0.1%. By inspection of the simulation 
data set as a whole, an even better agreement with experiment can be assumed. 
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5 Conclusion 



Transport and mixture properties for methanol and ethanol were predicted by 
molecular dynamics simulation on the basis of rigid, non-polarizable molecular 
models from preceding work. Self- and Maxwell- Stefan diffusion coefficients as 
well as shear viscosity were calculated by EMD and the Green-Kubo formalism. 
Comparing to experimental data for the two pure alcohols, present self-diffusion 
coefficient and shear viscosity data are very good to excellent over a temperature 
range of about 100 K. Slightly higher deviations were found for the shear vis- 
cosity at low temperatures, which are mainly due to the long-time behavior of 
the shear viscosity autocorrelation function. Regarding the mixture methanol -l- 
ethanol at ambient conditions over the full composition range, both self -diffusion 
coefficients also excellently agree with experimental data. The same holds for 
shear viscosity, especially in the methanol-rich region. Maxwell-Stefan diffusion 
coefficient data is, to our knowledge, not available in the literature. Present pre- 
dictions correspond very well with the model of Darken. The thermal conductivity 
was determined by BD reverse NEMD. Over a temperature range of about 100 K, 
an excellent agreement with the experiment was found as well for both pure fluids. 
Furthermore, predicted VLE data for the mixture methanol -I- ethanol along two 
isotherms correspond with experimental data practically throughout within the 
simulation uncertainty. The same holds for excess volume and excess enthalpy at 
ambient conditions in the full composition range. It has been demonstrated that 
present rigid, non-polarizable models for ethanol and methanol show a better per- 
formance to reproduce transport properties when compared to standard rigid and 
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flexible models available in the literature. 
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Appendix Simulation details 



In this work, EMD simulations for transport properties were made in two steps. 
In the first step, a short simulation in the isobaric-isothermal (NpT) ensemble was 
performed at ambient pressure to calculate the density at the desired temperature. 
In the second step, a canonic (NVT) ensemble simulation was performed at this 
temperature and density, to determine the transport properties. The simulations 
were carried out in a cubic box with periodic boundary conditions containing 
2048 molecules. Newton's equations of motion were solved using a fifth-order 
Gear predictor-corrector numerical integrator. The temperature was controlled 
by velocity scaling. In all EMD simulations, the integration time step was 0.98 
fs. The cut-off radius was set to rc = 15 A. Electrostatic long-range corrections 
were made using the reaction field technique with conducting boundary conditions 
{£Rf = oo) . On the basis of a center of mass cut-off scheme, the LJ long-range 
interactions were corrected using angle averaging [|63l . The simulations were 
equilibrated in the NVT ensemble over 10^ time steps, followed by production 
runs of 0.5 to 2 x 10^ time steps. Diffusion coefficients and shear viscosity were 
calculated using Eqs. (0]) to ([7]) with up to 9 900 independent time origins of the 
autocorrelation functions. The sampling length of the autocorrelation functions 
varied between 9 and 14 ps. The separation between time origins was chosen so 
that all autocorrelation functions have decayed at least to l/e of their normalized 
value to guarantee their time independence [64] . The uncertainties of the predicted 
values were estimated using a block averaging method (651. 

In the NEMD simulations with the PeX algorithm for the thermal conductiv- 
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ity, 800 molecules were placed in a parallelepiped box with = 3 = 3 Ly, 
where L, is the length of the simulation box in the different spatial directions. Pe- 
riodic boundary conditions were applied in all directions. The system was then 
equilibrated over 1 ns at the desired temperature and pressure by NpT simulation 
using a weak coupling bath [|66ll with long-range corrections [[65ll for pressure and 
energy. The resulting density of the equilibration run was then employed to gen- 
erate a new set of simulations in order to develop the thermal gradient in a run of 
2 ns using the NEMD scheme. In this case an exchange frequency X) = 3.3 fs^^ 
(i.e. every 150 time steps) for the energy flux was used. Electrostatic interactions 
were treated with the reaction field with conducting boundary conditions. The 
thermal conductivity was then computed as the average of 4 different runs over 1 
ns. The integration of the equations of motion was performed with a time step of 
2 fs with a cut-off radius of 12 A. A Verlet neighbor list was employed to improve 
the performance of the simulations. 

The VLB calculations of the mixture methanol -i- ethanol were performed us- 
ing the Grand Equilibrium method ll67ll in combination with the gradual insertion 
technique Il68l . The latter, a Monte Carlo based expanded ensemble method, was 
applied to determine the chemical potential in the liquid phase. In comparison to 
Widom's test particle method ll69l . where a full molecule is inserted into the fluid, 
the gradual insertion method introduces a fluctuating molecule. This molecule un- 
dergoes changes in a set of discrete states of coupling with all other real molecules 
of the fluid. Preferential sampling was performed in the vicinity of the fluctuating 
molecule to improve accuracy. The gradual insertion simulations were done with 
500 molecules in the liquid phase. Starting from a face-centered cubic lattice ar- 
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rangement, the simulation runs were equilibrated over 2 000 cycles in the NpT 
ensemble. The production phase was performed over 30 000 Monte Carlo cycles. 
Further simulation parameters of these runs were taken from Vrabec et al. f,70l . 

To evaluate the excess properties, EMD simulations in the NpT ensemble 
were performed. The system contained 1372 molecules in a cubic box with pe- 
riodic boundary conditions. The molecular trajectories were sampled with the 
algorithms and parameters as described above. The simulations were equilibrated 
in the NpT ensemble over 80 000 time steps, followed by a production run over 
600 000 time steps. 
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Table 1 : Literature review regarding transport properties of methanol, ethanol and 
their binary mixture calculated by molecular simulation. 





Self-diffusion 


Shear 
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coefficient 
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conductivity 
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Ethanol 


|[D|251|27l|2i|29l[30l|47l| 


mm 


mm 


Methanol + Ethanol 
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Table 2: Lennard- Jones, point charge and geometric parameters of the molecular 
models for methanol and ethanol, cf. Eq. ([U). Boltzmanns's constant is and the 
electronic charge is e. 



Site 










A 


K 


e 


Methanol 


ScH3 


3.7543 


120.592 


+0.24746 


SoH 


3.0300 


87.879 


-0.67874 


Sh 






+0.43128 


Ethanol 


ScH3 


3.6072 


120.15 




ScH2 


3.4612 


86.291 


+0.25560 


SoH 


3.1496 


85.053 


-0.69711 


Sh 






+0.44151 
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Table 3: Density, self-diffusion coefficient, shear viscosity and theraial conduc- 
tivity of pure liquid methanol and ethanol at 0. 1 MPa from present molecular dy- 
namics simulation. The number in parenthesis indicates the statistical uncertainty 
in the last digit. 



T 


P 










X 


K 


mol 


10-^m- 


-2s-l 


10-3pa s 


Wm-iR-i 


Methanol 


213 


27.02 


0.318 


(5) 


3.1 


(4) 




220 


26.80 


0.433 


(6) 


2.0 


(2) 


0.22 (3) 


240 


26.22 


0.750 


(8) 


1.7 


(2) 




260 


25.63 


1.21 


(1) 


1.0 


(2) 


0.207 (9) 


278.15 


25.14 


1.72 


(1) 


0.64 


(8) 




288 


24.81 


2.09 


(2) 


0.62 


(7) 




298.15 


24.51 


2.41 


(2) 


0.57 


(6) 




300 


24.50 










0.191 (9) 


308.15 


24.19 


2.97 


(2) 


0.47 


(5) 




318.15 


23.91 


3.54 


(2) 


0.41 


(6) 




320 


23.88 










0.188 (8) 


328.15 


23.58 


4.15 


(3) 


0.42 


(5) 




340.15 


23.21 


4.91 


(3) 


0.29 


(5) 




Ethanol 


239 


18.11 


0.192 


(6) 


3.1 


(3) 




253.15 


17.99 


0.250 


(5) 


2.4 


(2) 


0.18 (1) 


263 


17.83 


0.363 


(7) 


2.2 


(3) 




273.15 


17.62 


0.493 


(8) 


1.7 


(2) 


0.178(9) 


280 


17.46 


0.64 


(1) 


1.4 


(2) 




298.15 


17.08 


1.07 


(1) 


1.1 


(1) 


0.171 (9) 


308.15 


16.88 


1.33 


(2) 


1.0 


(2) 




320 


16.66 


1.68 


(1) 


0.8 


(1) 




328.15 


16.45 


2.01 


(2) 


0.61 


(8) 


0.161 (8) 


333 


16.31 


2.20 


(2) 


0.60 


(7) 




338.15 


16.26 


2.40 


(2) 


0.47 


(7) 
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Table 4: Density, self-diffusion and Maxwell-Stefan diffusion coefficient as well 
as shear viscosity of the mixture methanol (1) -i- ethanol (2) at 298.15 K and 0.1 
MPa from present molecular dynamics simulation. The number in parenthesis 
indicates the statistical uncertainty in the last digit. 





P 




D2 


£>12 




mol mol~^ 


mol 1-1 


lO-'^m-^s-i 


lO-'^m-^s-i 


lO-^m-V^ 


10-3pa s 





27.02 




1.07 (1) 






0.143 


26.22 


1.37 (2) 


1.12(1) 


1.4 (3) 


1.17(9) 


0.305 


25.63 


1.52 (1) 


1.25(1) 


1.8 (3) 


0.76 (9) 


0.500 


25.14 


1.73 (1) 


1.44(1) 


1.7 (2) 


0.62 (8) 


0.570 


24.81 


1.82 (1) 


1.53(1) 


1.6 (3) 


0.64 (9) 


0.754 


24.51 


2.07 (1) 


1.72 (2) 


2.3 (3) 


0.70 (9) 


0.796 


23.91 


2.14 (2) 


1.82 (2) 


2.0 (3) 


0.59 (7) 


0.890 


23.58 


2.28 (1) 


1.90 (2) 


2.6 (3) 


0.60 (8) 


1 


23.21 


4.95 (3) 
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Table 5: Vapor-liquid equilibrium data for the mixture methanol (1) + ethanol 
(2) at 393.15 K and 433.15 K from present molecular dynamics simulation. The 
number in parenthesis indicates the statistical uncertainty in the last digit. 



XI 


P 




yi 




P' 


P" 




mol mol^^ 


MPa 


mol mol 


mol 1-1 


mol 1-1 


kJ mol-i 


r = 393.15 K 





0.551 (6) 







14.84 (2) 


0.200(1) 


33.24 (7) 


0.150 


0.448 (9) 


0.207 


(5) 


15.62 (2) 


0.154 (3) 


33.86 (6) 


0.282 


0.482 (8) 


0.369 


(2) 


16.35 (2) 


0.168 (3) 


33.32 (7) 


0.492 


0.538 (9) 


0.589 


(2) 


17.54 (3) 


0.185 (3) 


32.67 (8) 


0.678 


0.588 (9) 


0.752 


(5) 


18.87 (3) 


0.212(3) 


31.92 (9) 


0.838 


0.572 (9) 


0.884 


(3) 


19.94 (3) 


0.205 (3) 


31.43 (8) 


1.0 


0.767 (5) 


1.0 




21.30 (3) 


0.294 (1) 


29.74 (8) 


r = 433.15 K 





1.32 


(3) 







13.58 (4) 


0.469 (2) 


28.41 (7) 


0.150 


1.29 


(2) 


0.192 


(4) 


14.34 (4) 


0.453 (7) 


28.52 (9) 


0.282 


1.38 


(2) 
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19.42 (4) 


0.711 (2) 


26.12(9) 
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Table 6: Volume, excess volume, enthalpy and excess enthalpy of the mixture 
methanol (1) + ethanol (2) at 298.15 K and 0.1 MPa from present molecular dy- 
namics simulation. The numbers in parenthesis indicate the statistical uncertainty 
in the last digits. 



xi V h 



ol mol ^ 


—3 

cm 


mol ^ 


_3 

cm 


mol ^ 


kJ mol-i 


J mol~^ 





58.44 


(5) 






-43.46 


(4) 




0.070 


57.15 


(5) 


0.00 


(7) 


-43.20 


(4) 


3 (52) 


0.200 


54.94 


(5) 


0.01 


(6) 


-42.75 


(4) 


2 (49) 


0.270 


53.72 


(5) 


0.03 


(6) 


-42.50 


(4) 


-1 (47) 


0.350 
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0.00 
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5 (44) 
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Abstract 

Transport properties of liquid methanol and ethanol are predicted by mole- 
cular dynamics simulation. The molecular models for the alcohols are rigid, 
non-polarizable and of united-atom type. They were developed in preced- 
ing work using experimental vapor-liquid equilibrium data only. Self- and 
Maxwell-Stefan diffusion coefficients as well as the shear viscosity of metha- 
nol, ethanol and their binary mixture are determined using equilibrium molec- 
ular dynamics and the Green-Kubo formalism. Non-equiUbrium molecular 
dynamics is used for predicting the thermal conductivity of the two pure 
substances. The transport properties of the fluids are calculated over a wide 
temperature range at ambient pressure and compared with experimental and 
simulation data from the literature. Overall, a very good agreement with the 
experiment is found. For instance, the self-diffusion coefficient and the shear 
viscosity are predicted with average deviations of less 8% for the pure alco- 
hols and 12% for the mixture. The predicted thermal conductivity agrees 
in average within 5% with the experimental data. Additionally, some ve- 
locity and shear viscosity autocorrelation functions are presented and dis- 
cussed. Radial distribution functions for ethanol are also presented. The 
predicted excess volume, excess enthalpy and the vapor-liquid equilibrium 
of the binary mixture methanol + ethanol are assessed and the vapor-liquid 
equiUbrium agree well with experimental data. 

Keywords: Green-Kubo; reverse NEMD; self-diffusion coefficient; Maxwell- 
Stefan diffusion coefficient; shear viscosity; thermal conductivity. 
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1 Introduction 



There has been significant progress in recent years to study and predict both static 
and dynamic thermophysical properties of fluids by molecular simulation. This is 
particularly interesting for transport properties of liquids, where, due to the com- 
plexity of the involved physical mechanisms, other theoretical approaches are of- 
ten unsatisfactory. Hence, molecular dynamics is a promising method for predict- 
ing transport properties of liquids. However, the number of publications regarding 
transport properties for technically relevant liquid systems that achieve quantita- 
tively accurate results is still low. The majority of these publications deals with the 
self-diffusion coefficient, leaving aside the more demanding Maxwell-Stefan dif- 
fusion coefficient, shear viscosity and thermal conductivity. Furthermore, molec- 
ular simulation studies on transport properties for complex liquids, e.g. mixtures 
containing highly polar hydrogen-bonding molecules, are rare. This is also the 
case for methanol, ethanol and their binary mixture. It should be pointed out 
that there is a series of articles on transport properties of pure liquid ethanol by 
Petravic et al. [?, ?, ?]. There are several works on transport properties of associ- 
ating mixtures: Typically, aqueous solutions of methanol [?, ?, ?, ?, ?] or ethanol 
[?, ?, ?, ?, ?] have been studied using molecular dynamics simulation. However, 
to our knowledge the mixture methanol + ethanol has not yet been regarded in that 
sense. 

Two fundamentally different methods for calculating transport properties by 
molecular dynamics simulation are available. The equilibrium methods (EMD), 
using either the Green-Kubo formalism or the Einstein relations, analyze the time 
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dependent response of a fluid system to spontaneous fluctuations. With non- 
equilibrium molecular dynamics (NEMD), on the other hand, the system's re- 
sponse to an externally applied perturbation is analyzed. NEMD was developed 
in order to increase the signal to noise ratio and to improve statistics and con- 
vergence. Both methods have different advantages and disadvantages, but are 
comparable in efficiency, as shown e.g. by Dysthe et al. [?]. 

The success of molecular dynamics simulation to predict thermodynamic prop- 
erties is highly determined by the potential model used to describe the intermolec- 
ular interactions. Numerous molecular models for methanol and ethanol have 
been proposed in the literature [?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?]. However, 
the ability of these models for predicting transport properties has been probed 
predominantly regarding the self -diffusion coefficient, cf. Table ??. 

For methanol, mainly the rigid OPLS molecular model from Jorgensen [?] and 
its modifications proposed by Haughney et al. [?] as well as by van Leeuwen and 
Smit [?], have been used to predict the self-diffusion coefficient. It was shown 
by comparison of the models by Jorgensen [?] and by Haughney et al. [?] that 
the rigid OPLS model and its modifications perform well for the self -diffusion 
coefficient, deviations to experimental data range from 7 to 20% [?]. Van de 
Ven-Lucassen et al. [?] showed that the model from van Leeuwen and Smit [?] 
predicts the self-diffusion coefficient of methanol at ambient temperature more 
accurately than the polarizable model of Caldwell and KoUman [?], and than the 
models by Haughney et al. [?]. Wheeler et al. [?] predicted the shear viscosity 
of methanol within 10% of experimental values. The thermal conductivity of 
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methanol at ambient conditions was reported by Nieto-Draghi [?] with a deviation 
of +11%. 

Different rigid [?, ?], semi-flexible [?], flexible [?] and polarizable [?, ?] 
molecular models for ethanol have been used to predict the self-diffusion coef- 
ficient. Thereby, it has been found that rigid, united-atom molecular models pre- 
dict the self-diffusion coefficient more reliably than polarizable, flexible models. 
Most ethanol models assessed by Wang and Cann [?] as well as the Jorgensen 
model [?] investigated by Saiz et al. [?] overestimate the self-diffusion coefficient 
rather significantly with deviations ranging from +8 to +50%. Allowing a po- 
tential model for bending results in a further increase in the predicted values [?]. 
Gonzalez et al. [?] found that rigid, non-polarizable models perform better than 
the polarizable ones for predictions of the self -diffusion coefficient. The shear vis- 
cosity of ethanol has been predicted using EMD [?] and NEMD [?]. Both works 
achieved good results (average deviations are 11 and 7%, respectively) using the 
united-atom model and the all-atom model of Jorgensen et al. [?, ?], respectively. 
Simulations performed by Petravic [?] underestimate the thermal conductivity of 
liquid ethanol near ambient temperatures by approximately —10%. 

The present article assesses the predictive power of rigid, non-polarizable, 
united-atom type Lenard- Jones based models with superimposed point charges 
for methanol and ethanol developed in prior work of our group [?, ?]. These 
models have been optimized on the basis of experimental data of vapor pressure 
and saturated liquid density in the whole temperature range from the triple point 
to the critical point. No experimental data for transport properties was thereby 
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taken into account. The most significant transport properties, i.e. self-diffusion 
coefficient, shear viscosity and thermal conductivity, were regarded here for the 
pure liquids using EMD and NEMD. Furthermore, the two molecular models were 
tested regarding the prediction of self- and Maxwell-Stefan diffusion coefficients 
and shear viscosity of the binary mixture. Additionally, binary vapor-liquid equi- 
libria as well as excess volume and excess enthalpy were predicted. The goal of 
this study is to demonstrate the ability of rigid, non-polarizable molecular models, 
adjusted to pure substance vapor-liquid equilibria only,to accurately describe the 
transport properties of strongly polar and H-bonding liquids. 

This work is organized as follows: Firstly, the employed molecular models 
are briefly described. Then, the simulation techniques are explained, the results 
for self -diffusion coefficient, shear viscosity and thermal conductivity of the pure 
liquids are presented and compared to experimental and simulation data from the 
literature. The self- and Maxwell-Stefan diffusion coefficients as well as the shear 
viscosity for the mixture methanol -I- ethanol are presented and compared with 
experimental data where possible. Furthermore, the velocity and shear viscosity 
autocorrelation functions are described and discussed. Simulated radial distribu- 
tion functions of ethanol are also presented. Thereafter, the simulation results for 
vapor-liquid equilibria, excess volume and excess enthalpy for the binary mixture 
are given in comparison to experimental values. Finally, conclusions are drawn. 
Simulation details are summarized in the Appendix. 
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2 Molecular models 



Throughout this work, rigid, non-polarizable molecular models of united-atom 
type from earlier work of our group [?, ?] were used. Both models account for the 
intermolecular interactions, including hydrogen bonding, by a set of united-atom 
Lennard- Jones (LJ) sites and superimposed point charges. Hence, the potential 
energy Uij between two alcohol molecules / and j can be written as 

n m 

Uij {njab) = 12 12 "^^"^ 
a=lb=l 

where a is the site index of molecule /, b the site index of molecule j, and n 
and m are the numbers of interaction sites of molecules i and j, respectively, ryab 
represents the site-site distance between molecules / and j. The LJ size and energy 
parameters are Oab and Eab- lia and qjb are the point charges located at the sites a 
and b on the molecules / and j, and 8o is the permittivity of vacuum. 

The methanol model has two LJ sites, i.e. one for the methyl group and one 
for the hydroxyl group. Additionally, three point charges are superimposed, i.e. 
one at each of the LJ site centers and one at the nucleus position of the hydroxyl 
hydrogen. Ethanol was modeled by three LJ sites, i.e. one for the methyl, the 
methylene and the hydroxyl group each. Three point charges are located here as 
well on the methylene and hydroxyl LJ site centers and on the nucleus position of 
the hydroxyl hydrogen. 

The parameter set of the methanol model was obtained using the optimization 
procedure of StoU [?] on the basis of the geometry of the model from van Leeuwen 



^ijabj \fijabj 



+ 



QiaQ jb 
4n£orijab ' 



(1) 
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and Smit [?]. The geometry of the ethanol model was derived from ab initio 
quantum chemistry calculations. The AUA4 parameters of Ungerer et al. [?] 
were applied for the methyl and methylene LJ sites. The point charge, LJ size 
and energy parameters as well as the offset of the hydroxyl LJ site were fitted 
to experimental vapor pressure and saturated liquid density. It should be pointed 
out that no information on transport or mixture properties was included in the 
optimization procedures of any molecular model employed here so that the data 
presented below is strictly predictive. The model parameters were taken from 
[?, ?] and are summarized in Table ??. 

To define a molecular model for a binary mixture on the basis of pairwise ad- 
ditive pure substance models, only the unlike interactions have to be specified. In 
case of polar interaction sites, i.e. point charges here, this can straightforwardly be 
done using the laws of electrostatics. However, for the unlike LJ parameters there 
is no physically sound approach [?], and for predictions combination rules have 
to be employed. Here, the interaction between unlike LJ sites of two molecules 
are calculated applying the Lorentz-Berthelot combining rules: 



2 



(2) 



and 



^ab — \/ ^aa^bb ■ 



(3) 
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3 Methodology 



Dynamic properties can be obtained from EMD simulations by means of the 
Green-Kubo formalism [?, ?]. These equations are a direct relationship between a 
transport coefficient and the time integral of an autocorrelation function of a par- 
ticular microscopic flux in a system in equilibrium. This method was used in this 
work to calculate self- and Maxwell-Stefan diffusion coefficients as well as the 
shear viscosity. The thermal conductivity was calculated using a modified version 
of the boundary driven (BD) reverse NEMD algorithm from Miiller-Plathe [?], 
also referred to as PeX algorithm. 

3.1 Diffusion coefficients 

The self-diffusion coefficient Di is related to the mass flux of single molecules 
within a fluid. Therefore, the relevant Green-Kubo expression is based on the 

individual molecule velocity autocorrelation function as follows 



where Vfc(?) is the center of mass velocity vector of molecule k at some time t, and 
<...> denotes the ensemble average. Eq. (??) is an average over all molecules 
in a simulation, since all contribute to the self -diffusion coefficient. In a binary 
mixture, the Maxwell-Stefan diffusion coefficient Dij can also be written in terms 
of the center of mass velocity 




(4) 
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where Mi and Xj are the molar mass and mole fraction of species i. From the 
expression above, the collective character of the Maxwell-Stefan diffusion coeffi- 
cient is evident. 

Since the present simulations provide self- as well as the Maxwell-Stefan dif- 
fusion coefficients simultaneously, a comparison to the simple predictive approach 
suggested by Darken [?] is possible. The Maxwell-Stefan diffusion coefficient D,^ 
is estimated from the self -diffusion coefficients of the two components Z),- and Dj 
in the binary mixture [?] 

Dij^Xi-Di + XjDj- (6) 

This approach was found to be superior to the correlations by Vignes [?] and 
by Caldwell and Babb [?] in prior work on three binary mixtures [?]. However, it 
is of little practical use due to the fact that it relies on the self-diffusion coefficients 
in the mixture, which are rarely available. 

3.2 Shear viscosity 

The shear viscosity r|, as defined by Newton's "law" of viscosity, is a measure of 
the resistance of a fluid to a shearing force [?]. It is associated with momentum 
transport under the influence of velocity gradients. Hence, the shear viscosity can 
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be related to the time autocorrelation function of the off-diagonal elements of the 
stress tensor J^p [?] 

"^^Vk^L ^^(•^p'W-^P^W), (7) 

where V stands for the molar volume, ka is the Boltzmann constant, and T denotes 
the temperature. Averaging over all three independent elements of the stress ten- 
sor, i.e. , Jp^,andJp^, improves the statistics of the simulation. The component 
Jp' of the microscopic stress tensor Jp is given in terms of the molecule positions 
and velocities by [?] 



N -t N N n n 

^f=i'»vfvj-ii:EEE4^^. (8) 

!=1 ^ ;=iy^i/t=u=i ^f^kl 

Here, the lower indices / and k count the interaction sites, and the upper indices 
X and y denote the spatial vector components, e.g. for velocity v- or site-site dis- 
tance r-j. The first term of Eq. (??) is the kinetic energy contribution and the 
second term is the potential energy contribution to the shear viscosity. Conse- 
quently, the Green-Kubo integral (??) can be decomposed into three parts, i.e. 
the kinetic energy contribution, the potential energy contribution and the mixed 
kinetic-potential energy contribution [?]. Eqs. (??) and (??) may directly be ap- 
plied to mixtures. 
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3.3 Thermal conductivity 

The thermal conductivity X, as defined by Fourier's "law" of heat conduction, 
characterizes the capability of a substance for molecular energy transport driven 
by temperature gradients. In BD-NEMD simulation used in this work, two slabs of 
a given thickness, large enough to contain a sufficiently large number of molecules 
on average but much smaller than the total length of the simulation volume in z 
direction, are defined. In order to create a heat flux in z direction, the molecule 
with the highest kinetic energy in the cold slab is selected with a frequency D 
to perform a hypothetical elastic collision with the molecule having the lowest 
kinetic energy in the hot slab [?]. After such a virtual collision the two molecules 
have exchanged their momentum without changing their respective position in 
space. The new velocity of the molecule in the cold region is then 



old I ,„ ..old 



new_ „old , m,Vr+mhVl 

— "I \ ; W 

nic + mh 



and, for the molecule in the hot region 



^„^_^.,,^r^l±!!^ (10) 

rric + mh 

where and m/j are the respective masses of the selected molecules in the cold 
and hot slabs, v^''^ and stand, respectively, for the velocities of the molecules 
before and after the collision. 

During the elastic collision the molecules exchange their momentum. Thus, 

12 



the total momentum and the total energy of the system remain constant. The ki- 
netic energy exchanged by means of these elastic collisions determines the energy 
flux in the system. The heat flow density in the steady state is then given by 

^•^^^^=21^ ^ ^((vr)'-(vrf)- (11) 

transfers 

In this expression, A is the cross sectional area in x and y direction, t is the time 
interval of measure in the steady state during which a given number of momentum 
transfers have occurred. {Jc)( is the heat flux in z direction, which can indistinctly 
be evaluated in the cold or the hot slab. It is required that the number of molecules 
in both slabs is sufficiently large so that the intermolecular interactions within the 
slabs will properly thermalize the velocity distribution before a new collision takes 
place. The thermal gradient Vr^. is obtained from a linear fit of the temperature 
profiles from the simulation. In the steady state, the thermal conductivity can be 
obtained by 



4 Simulation results 
4.1 Diffusion coefficients 

The self -diffusion coefficient of pure methanol and pure ethanol was predicted at 
ambient pressure in the temperature range from 210 to 340 K. Present numerical 
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data are given in Table ??. Figures ?? and ?? show the self-diffusion coefficient 
for methanol and ethanol as a function of temperature, respectively. The statistical 
error of the simulation data is estimated to be in the order of 1%. 

The self-diffusion coefficient of methanol shows a very good agreement with 
the experimental values. At temperatures below 280 K, the self-diffusion coef- 
ficient was underpredicted by approximately —7%. Above ambient temperature, 
the mobility of the model molecules increases slightly more rapidly than that of 
real methanol molecules. Hence, for temperatures above 300 K the self-diffusion 
coefficient was overestimated on average by +8%. However, the significant scat- 
ter of the experimental data should be taken into account. In comparison to pre- 
vious simulation studies from the literature, the present results cover the widest 
temperature range. The broadest available study by Haughney et al. [?], using a 
modified Jorgensen model, is comparable with present results near ambient tem- 
perature where average deviations are 9 %. However, the data from Haughney et 
al. is scattering and the self -diffusion coefficient is significantly underestimated 
above 300 K. 

There is an excellent agreement between experimental data and present self- 
diffusion coefficient for ethanol. The predicted data also correctly reproduces 
the temperature dependence over the whole temperature range from 239 to 338 K. 
Contrary to previous simulation reports [?, ?, ?, ?, ?, ?, ?], where the self-diffusion 
coefficient of ethanol has mostly been significantly overestimated, deviations of 
the present data are negative and on average only —6%. This is explained by the 
different potential model used in this work. 
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Self -diffusion coefficients of methanol and ethanol in their binary mixture 
were predicted at ambient conditions for the entire composition range. Numer- 
ical results are presented in Table ?? and plotted in Figure ??. The estimated 
statistical uncertainty of the simulation data is also 1%. As shown in Figure ??, 
there is an excellent agreement between predicted self-diffusion coefficients and 
the experimental data for both alcohols. Present simulation results slightly over- 
estimate the self-diffusion of both alcohols by +5% on average. Furthermore, the 
predicted data reproduce the composition dependence well. 

The Maxwell-Stefan diffusion coefficient of the mixture methanol + ethanol 
calculated by molecular simulation is reported in Table ??. In Figure ??, the pre- 
dicted values are plotted as a function of the composition. The statistical uncer- 
tainties are on average 15%. These large error bars are attributed to the collective 
nature of the Maxwell-Stefan diffusion coefficient. Due to the lack of experi- 
mental data, present simulation results are only compared to the predictions of 
Darken's model [?], cf. Eq. (??). This model can be considered to be a good 
approximation here, in agreement to previous findings for three other mixtures 
[?]. 

4.2 Shear viscosity 

The shear viscosity of pure methanol and pure ethanol was predicted at ambient 
pressure for a temperature range of about 100 K around ambient conditions and is 
reported in Table ??. The temperature dependence of the shear viscosity is shown 
in Figures ?? and ?? for methanol and ethanol, respectively. Statistical errors are 
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in the range from 10 to 17%. 

In the case of methanol, the predicted shear viscosity agrees very well with 
the experimental data over the whole regarded temperature range with an average 
deviation of 8%. At higher temperatures, the predicted viscosities are closer to the 
experimental ones than at lower temperatures, cf. Figure ??. It should be pointed 
out that the highest deviation, being 20%, was found at the lowest regarded tem- 
perature 220 K. Here, the high density in combination with strongly interacting 
molecules causes large simulation uncertainties due to the long-time behavior of 
the shear viscosity autocorrelation function. Present shear viscosity predictions 
are better than the EMD data of Wheeler et al. [?] and are comparable with the 
NEMD Edwald data of the same author. 

The calculated shear viscosity for ethanol shows a good agreement to exper- 
imental data with an average deviation of 8%. It can be noticed in Figure ?? 
that there is a tendency to underestimate the shear viscosity. Petravic and Del- 
hommelle [?] also underestimated the shear viscosity using the OPLS model from 
Jorgensen [?]. At high densities, however, present data better predict the shear 
viscosity. Zhao et al. [?] recently reported shear viscosity data using the OPLS- 
AA model [?] and NEMD which is in very good agreement with experimental 
data. 

The three contributions to the shear viscosity were also calculated. In agree- 
ment with the results of Meier et al. [?], the potential energy contribution is domi- 
nant for both studied alcohols, r\pp > 90%. In general, for the studied temperature 
range, the kinetic energy term has a negligible positive contribution to the shear 
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viscosity, 'r\kk < 2%. The mixed kinetic-potential energy term may contribute pos- 
itively or negatively to the shear viscosity. Here, this contribution varies between 
0.1 and 8%. 

The simulated shear viscosity data of the mixture methanol + ethanol is given 
in Table ?? and plotted in Figure ?? together with the available experimental data. 
There is a very good agreement between experimental and predicted shear vis- 
cosity, being mostly within the simulation uncertainty of around 12%. Deviations 
from the experimental data are on average 12% as well. 

4.3 Thermal conductivity 

Present results for the thermal conductivity of pure methanol and pure ethanol at 
ambient pressure are summarized in Table ??. The agreement with the experi- 
mental data is remarkable with an average deviation of 5% for methanol and of 
2% for ethanol, cf. Figures ?? and ??. The thermal conductivity for methanol is 
slightly underestimated by the present predictions, while for ethanol it is slightly 
overestimated. Previous results for the thermal conductivity of methanol [?] and 
ethanol [?], using the OPLS models, show deviations of about 11 and 10%, re- 
spectively. In contrast to present results, the calculated thermal conductivity by 
Petravic [?] lies below the experimental data. 
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4.4 Time dependent behavior of the Green-Kubo integrands 
Velocity autocorrelation function 

The velocity autocorrelation function provides an interesting insight into the liq- 
uid state. Selected normalized velocity autocorrelation functions (VACF) of pure 
methanol and ethanol are shown in Figure ??. In agreement with findings in the 
literature [?], VACF may have two minima at short times (t < 0.5 ps). At low 
temperatures and high densities, the VACF of both pure liquids decrease sharply 
and become negative. After this first minimum, they increase slightly to drop 
again into a lower minimum, for methanol, or to a higher minimum, for ethanol. 
Beyond the second minimum, the VACF converge to zero following a character- 
istic path. The observed negative values of the VACF are related to backscat- 
tering collisions, also known as cage-effect, which govern short time dynamics. 
Here, the free movement of the molecules is hindered by surrounding molecules, 
which form a cage-like structure. Michels and Trappeniers [?] attributed this phe- 
nomenon to the formation of "bound states" since this effect can only be observed 
for molecular models that include attractive forces. 

For higher temperatures at ambient pressure, where the liquid density is lower, 
molecular collisions are less frequent and the depth of the minima is gradually 
reduced until they almost disappear, cf. Figure ??. At short times, the VACF decay 
faster at lower temperatures. Moreover, the initial decay of the VACF of methanol 
occurs earlier than in the corresponding VACF of ethanol. Methanol VACF exhibit 
a more pronounced backscattering oscillations than those of ethanol, particularly 
at low temperatures near its melting temperature (Tm = 175.37 K) [?]. 
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Figure ?? shows the behavior of the integrals of the VACF given by Eq. (??). 
It can be seen that the self-diffusion coefficient converges in all regarded cases to 
its final value after less than 4 ps. Thus, after this time interval the long-time tail 
of the VACF does not contribute significantly anymore. 

The behavior of the VACF of methanol and ethanol was also considered in the 
mixture. From Figures ?? and ?? it can be seen that the VACF of methanol and 
ethanol in the mixture differ from those of the pure alcohols at the same tempera- 
ture and pressure. Methanol VACF exhibit a significantly increased backscattering 
due to the presence of ethanol, cf. Figure ??. This indicates that at short times, 
ethanol molecules further restrain the mobility of methanol molecules. Ethanol 
VACF also show a composition dependent behavior. In this case, the backscat- 
tering is reduced as the concentration of methanol increases, cf. Figure ??. This 
suggests that methanol molecules prevent the formation of cage-like structures 
here, allowing ethanol molecules to move more freely in the mixture. 

Shear viscosity autocorrelation function 

Alder et al. [?] were the first to suggest the existence of a long-time tail in the 
shear viscosity autocorrelation function (SACF) near the phase boundary to solid- 
ification due to the low compressibility of such liquid states. Luo and Hoheisel 
[?, ?] showed that the long-time behavior of the SACF is determined primarily by 
long-lived orientational correlations, occurring in liquids containing anisotropic 
molecules. This behavior causes convergence problems in the Green-Kubo inte- 
gral, since the SACF must be calculated over a relatively long time interval. 
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Figure ?? shows selected normalized SACF of pure methanol. The presence 
of a significant long-time tail can clearly be noticed in Figure ?? where the inte- 
grals of the SACF of methanol are plotted for three temperatures. As expected, 
the behavior of the long-time tail of the SACF is highly dependent on molecule 
species and state point. For pure methanol at a low temperature of 260 K, approx- 
imately 60% of the contribution to the total shear viscosity comes from the first 
2 ps. At a higher temperature of 340 K the first 2 ps contribute 90% to the total 
shear viscosity. In the case of ethanol (not shown here), the SACF are significant 
for even longer times, i.e. at 263 K the first 2 ps of the SACF contribute only 45% 
to the final value of the shear viscosity. At 338 K this contribution reaches 75%. 
This difference can be explained by the higher molecular anisotropy of ethanol. 

The kinetic, potential and mixed kinetic-potential energy contributions to the 
SACF were also considered in this work. Over the whole regarded range of state 
points, the behavior of the kinetic contribution is analogous to the mono tonic 
chair-like decay of autocorrelation functions for non-attractive fluids. The cor- 
responding time integral exhibits no important contribution from the long-time 
behavior, cf. Figure ??. The potential energy contribution shows at short times 
oscillations similar to those of the VACF, cf. Figures ?? and ??. These oscillations 
have been also related to the presence of "bound states" by Michels and Trappe- 
niers [?] and by Meier et al. [?]. The long-time tail is significant and approximates 
the behavior of the total SACF. The mixed kinetic-potential energy contribution 
shows a completely different behavior at short times, cf. Figure ??, where data is 
presented for pure methanol at four temperatures. At very short times (? < 0. 1 ps), 
two roughly mirrored patterns can be seen, the SACF decay/raise sharply to reach 
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a minimum/maximum and oscillate strongly during their decay to zero. Meier et 
al. [?], who investigated simple LJ fluids, only observed initial positive peaks in 
the mixed kinetic-potential energy contribution. To our knowledge, the presence 
of negative initial peaks has not been reported before. The long-time tail behavior 
of the mixed contribution is difficult to analyze because of the presence of strong 
oscillations. 

4.5 Hydrogen bonding dynamics 

Petravic et al. [?] have shown that internal degrees of freedom play a crucial 
role on the relaxation of the hydrogen bonding dynamics of liquid ethanol. For 
instance, they found that "freezing" bending and torsion in the OPLS model of 
ethanol produces a slowdown of the dynamics of the hydrogen bonds at 273 K 
(they observed an increase of the hydrogen bond life time of about 27%). For the 
methanol model used in this work, an excellent agreement between simulated and 
experimental site-site radial distribution functions and hydrogen bond statistics 
has been reported in previous work [?]. In the present work, the hydrogen bond 
dynamics of the rigid ethanol model was compared with the results using a flexi- 
ble model version including bending and torsion. For this porpouse, a geometrical 
hydrogen bonding criterion with voh below 2.6 A, rgo below 3.5 A and Q(oh...o) 
below 30 ° together with a life time probability accounting for intermittent hy- 
drogen bonds [?] was used. In the flexible version, bending and torsion angles 
were allowed to evolve using a standard harmonic and a Fourier series potential, 
respectively. The potential constants were taken from the TraPPE intermolecu- 
lar potential for alcohols [?], keeping the equilibrium values of the bending and 
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torsion angles of the rigid model and all the remaining inter-molecular model pa- 
rameters as well as the bond distances between all sites. Simulations were then 
performed at a specified temperature of 273 K (using the Nose-Hoover thermo- 
stat) and a density of 17.62 mol/1 to compute the radial distribution functions, the 
hydrogen bond dynamics and self-diffusion coefficient for both, rigid and flexible 
version of the present ethanol model. 

The comparison of the radial distribution functions from the flexible and the 
rigid model can be seen in Figure ??. It was found that there are only minor differ- 
ences in the liquid structure as revealed by the three radial distribution functions 
goO', gOH and gHH- The position and height of the first peaks are almost identi- 
cal for both models, with the exception of gHH, where a more pronounced peak 
is observed for the rigid model, cf. Figure ??. Only small differences are noted 
around the second peak of goo and goH, mainly related to the molecules hydrogen 
bonded in the second solvation shell. In the latter case, a more pronounced second 
peak was found for the rigid model. In general, it can be said that the inclusion 
of the flexibility does not strongly affect the structure of liquid ethanol around the 
first solvation shell. 

In the next step, differences in the hydrogen bond dynamics for the two mod- 
els were analyzed. Average life times of hydrogen bonds of 137 and 80 ps were 
obtained for the rigid and flexible versions, respectively. The flexibility signifi- 
cantly reduces the hydrogen bond life time by a factor of 1.7. Petravic et al. found 
in a similar investigation for the OPLS ethanol model a reduction in the hydrogen 
bond dynamics of about 1.3. Hence, the effect is more pronounced for the present 
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model. The hydrogen bond dynamics has a significant impact on the transport 
properties. Eg. at 273 K a self-diffusion coefficient of 0.802 x 10~^m^/s was 
obtained for the flexible version, which is in perfect agreement with the values re- 
ported by Petravic et al. [?]. This value is higher by a factor of 1 .6 compared to the 
self -diffusion coefficient obtained for the rigid model (0.493 x 10~^m^/s) which 
is in excellent agreement with the experimental data (0.48 x 10~^m^/s) [?]. This 
comparison shows the impact of the internal degrees of freedom on the hydrogen 
bond dynamics and the self -diffusion coefficient. Similar effects can be expected 
also for other transport properties. Flexible and rigid models of H-Bonding sys- 
tems obviously show significant differences regarding transport properties. This 
should also be kept in mind when assessing the comparisons of different models 
shown throughout this work. 

4.6 Vapor-liquid equilibria 

Vapor-liquid equilibria (VLE) often serve as a benchmark to assess molecular 
models. This is a useful choice, due to the fact that VLE represent the dominant 
feature in the region of fluid states. Indeed, they are nowadays preferably being 
used in the development of molecular models by optimizing the model parameters 
to VLE properties. This approach was also followed in the development of both 
regarded pure fluid models. In the next step, molecular models can be assessed 
on how they perform in mixtures. Based on experience for numerous binaries and 
ternaries, very good success can be expected, if the pure substance models are 
good. 
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Calculated VLE data of pure methanol and ethanol using present potential 
models have been reported previously for a temperature range from 270 to 490 
K [?, ?]. The average deviations found for methanol in vapor pressure, saturated 
liquid density and heat of vaporization are of 1.1, 0.6, and 5.5%, respectively 
[?]. The corresponding reported average deviations for ethanol are of 3.7, 0.3 and 
0.9%, respectively [?]. 

Present simulation results for the binary VLE are given in Table ?? and com- 
pared in Figures ?? and ?? to experimental data by Butcher and Robinson [?] and 
to the Peng-Robinson equation of state [?]. As in the molecular simulations no 
adjustment to binary data was carried out, also the Peng-Robinson equation was 
used without any adjustment to such data (kij = 0). Figure ?? shows the almost 
"ideal" behavior of methanol + ethanol for two temperatures. In fact, these were 
the highest two isotherms for which experimental data were available to us. At the 
higher one (443. 15 K), however, it can be seen from the experimental data that the 
bubble line is not fully linear. The comparison of simulation and experiment is ex- 
cellent throughout, both the pressure slope and the composition relation between 
vapor and liquid match good, only minor deviations are found in the ethanol-rich 
composition range where slightly too high vapor compositions are obtained. 

4.7 Excess properties 

Excess properties are basic mixing properties which contain real mixing behavior. 
In case of methanol + ethanol non-ideal mixing effects are very weak. At equimo- 
lar composition, where the mixing influence is strongest, experimental data sug- 
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gest a volume increase of 0.017% only. Also for enthalpy, a tiny positive excess 
contribution of 0.011 % was found by experiment. This resolution is difficult to 
achieve by molecular simulation. 

To determine excess properties, a straightforward approach was used here. 
Three simulations at specified temperature and pressure were performed, two for 
the pure substances and one of the mixture at a given composition. A binary 
excess property given by 

^y-xi-yi-X2-y2, (13) 

where y^ represents the excess volume or the excess enthalpy h^. Enthalpies 
or volumes of the pure liquids methanol and ethanol as well as of their binary 
mixture are denoted by yi, y2 and y, respectively. 

Table ?? contains present simulation data for ambient conditions which are 
compared to experimental results [?, ?, ?, ?] in Figures ?? and ??. The excess 
volume, cf. Figure ??, agrees with the experiment in all cases within its seem- 
ingly large error bars. However, it should be noted that the simulation uncertainty 
is only around 0.1% of the total volume of the mixture. A slight tendency of the 
simulation data towards higher excess volumes might be assumed. In Figure ??, 
the present excess enthalpy is shown together with experimental data. The maxi- 
mal excess contribution from experiment is only 0.01 1%, which is well below the 
simulation uncertainty of around 0.1%. By inspection of the simulation data set 
as a whole, an even better agreement with experiment can be assumed. 
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5 Conclusion 



Transport and mixture properties for methanol and ethanol were predicted by 
molecular dynamics simulation on the basis of rigid, non-polarizable molecular 
models from preceding work. Self- and Maxwell-Stefan diffusion coefficients as 
well as shear viscosity were calculated by EMD and the Green-Kubo formalism. 
Comparing to experimental data for the two pure alcohols, present self-diffusion 
coefficient and shear viscosity data are very good to excellent over a temperature 
range of about 100 K. Slightly higher deviations were found for the shear vis- 
cosity at low temperatures, which are mainly due to the long-time behavior of 
the shear viscosity autocorrelation function. Regarding the mixture methanol + 
ethanol at ambient conditions over the full composition range, both self -diffusion 
coefficients also excellently agree with experimental data. The same holds for 
shear viscosity, especially in the methanol-rich region. Maxwell-Stefan diffusion 
coefficient data is, to our knowledge, not available in the literature. Present pre- 
dictions correspond very well with the model of Darken. The thermal conductivity 
was determined by BD reverse NEMD. Over a temperature range of about 100 K, 
an excellent agreement with the experiment was found as well for both pure fluids. 
Furthermore, predicted VLE data for the mixture methanol + ethanol along two 
isotherms correspond with experimental data practically throughout within the 
simulation uncertainty. The same holds for excess volume and excess enthalpy at 
ambient conditions in the full composition range. It has been demonstrated that 
present rigid, non-polarizable models for ethanol and methanol show a better per- 
formance to reproduce transport properties when compared to standard rigid and 
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flexible models available in the literature. 
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Appendix Simulation details 



In this work, EMD simulations for transport properties were made in two steps. 
In the first step, a short simulation in the isobaric-isothermal (NpT) ensemble was 
performed at ambient pressure to calculate the density at the desired temperature. 
In the second step, a canonic (NVT) ensemble simulation was performed at this 
temperature and density, to determine the transport properties. The simulations 
were carried out in a cubic box with periodic boundary conditions containing 
2048 molecules. Newton's equations of motion were solved using a fifth-order 
Gear predictor-corrector numerical integrator. The temperature was controlled 
by velocity scaling. In all EMD simulations, the integration time step was 0.98 
fs. The cut-off radius was set to = 15 A. Electrostatic long-range corrections 
were made using the reaction field technique with conducting boundary conditions 
{^RF = °°). On the basis of a center of mass cut-off scheme, the LJ long-range 
interactions were corrected using angle averaging [?]. The simulations were equi- 
librated in the NVT ensemble over 10^ time steps, followed by production runs 
of 0.5 to 2 X 10^ time steps. Diffusion coefficients and shear viscosity were cal- 
culated using Eqs. (??) to (??) with up to 9 900 independent time origins of the 
autocorrelation functions. The sampling length of the autocorrelation functions 
varied between 9 and 14 ps. The separation between time origins was chosen so 
that all autocorrelation functions have decayed at least to l/e of their normalized 
value to guarantee their time independence [?]. The uncertainties of the predicted 
values were estimated using a block averaging method [?]. 

In the NEMD simulations with the PeX algorithm for the thermal conductiv- 
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ity, 800 molecules were placed in a parallelepiped box with Lj. = 3 = 3 Ly, 
where L, is the length of the simulation box in the different spatial directions. Pe- 
riodic boundary conditions were applied in all directions. The system was then 
equilibrated over 1 ns at the desired temperature and pressure by NpT simulation 
using a weak coupling bath [?] with long-range corrections [?] for pressure and 
energy. The resulting density of the equilibration run was then employed to gen- 
erate a new set of simulations in order to develop the thermal gradient in a run of 
2 ns using the NEMD scheme. In this case an exchange frequency D = 3.3 fs~^ 
(i.e. every 150 time steps) for the energy flux was used. Electrostatic interactions 
were treated with the reaction field with conducting boundary conditions. The 
thermal conductivity was then computed as the average of 4 different runs over 1 
ns. The integration of the equations of motion was performed with a time step of 
2 fs with a cut-off radius of 12 A. A Verlet neighbor list was employed to improve 
the performance of the simulations. 

The VLE calculations of the mixture methanol -I- ethanol were performed us- 
ing the Grand Equilibrium method [?] in combination with the gradual insertion 
technique [?]. The latter, a Monte Carlo based expanded ensemble method, was 
applied to determine the chemical potential in the liquid phase. In comparison to 
Widom's test particle method [?], where a full molecule is inserted into the fluid, 
the gradual insertion method introduces a fluctuating molecule. This molecule un- 
dergoes changes in a set of discrete states of coupling with all other real molecules 
of the fluid. Preferential sampling was performed in the vicinity of the fluctuating 
molecule to improve accuracy. The gradual insertion simulations were done with 
500 molecules in the liquid phase. Starting from a face-centered cubic lattice ar- 
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rangement, the simulation runs were equilibrated over 2 000 cycles in the NpT 
ensemble. The production phase was performed over 30 000 Monte Carlo cycles. 
Further simulation parameters of these runs were taken from Vrabec et al. [?]. 

To evaluate the excess properties, EMD simulations in the NpT ensemble 
were performed. The system contained 1372 molecules in a cubic box with pe- 
riodic boundary conditions. The molecular trajectories were sampled with the 
algorithms and parameters as described above. The simulations were equilibrated 
in the NpT ensemble over 80 000 time steps, followed by a production run over 
600 000 time steps. 
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Table 1 : Literature review regarding transport properties of methanol, ethanol and 
their binary mixture calculated by molecular simulation. 





Self-diffusion Shear 


Thermal 




coefficient viscosity 


conductivity 


Methanol 


[7 ? ? 7 7J [7] 


[?] 


Ethanol 


[7777777] ['^ 


[?, ?] 


Methanol + Ethanol 
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Table 2: Lennard- Jones, point charge and geometric parameters of the molecular 
models for methanol and ethanol, cf. Eq. (??). Boltzmanns's constant is and 
the electronic charge is e. 



Site 










A 


K 


e 


Methanol 


SCH3 


3.7543 


120.592 


+0.24746 


SoH 


3.0300 


87.879 


-0.67874 


Sh 






+0.43128 


Ethanol 


ScH3 


3.6072 


120.15 




ScH2 


3.4612 


86.291 


+0.25560 


SoH 


3.1496 


85.053 


-0.69711 


Sh 






+0.44151 
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Table 3: Density, self-diffusion coefficient, shear viscosity and thermal conduc- 
tivity of pure liquid methanol and ethanol at 0. 1 MPa from present molecular dy- 
namics simulation. The number in parenthesis indicates the statistical uncertainty 
in the last digit. 



T 


P 








X 


K 


mol 1-1 


lO-^m" 


-2s-l 


10-3pa s 


Wm-iR-i 


Methanol 


213 


27.02 


0.318 


(5) 


3.1 (4) 




220 


26.80 


0.433 


(6) 


2.0 (2) 


0.22 (3) 


240 


26.22 


0.750 


(8) 


1.7 (2) 




260 


25.63 


1.21 


(1) 


1.0 (2) 


0.207 (9) 


278.15 


25.14 


1.72 


(1) 


0.64 (8) 




288 


24.81 


2.09 


(2) 


0.62 (7) 




298.15 


24.51 


2.41 


(2) 


0.57 (6) 




300 


24.50 








0.191 (9) 


308.15 


24.19 


2.97 


(2) 


0.47 (5) 




318.15 


23.91 


3.54 


(2) 


0.41 (6) 




320 


23.88 








0.188(8) 


328.15 


23.58 


4.15 


(3) 


0.42 (5) 




340.15 


23.21 


4.91 


(3) 


0.29 (5) 




Ethanol 


239 


18.11 


0.192 


(6) 


3.1 (3) 




253.15 


17.99 


0.250 


(5) 


2.4 (2) 


0.18 (1) 


263 


17.83 


0.363 


(7) 


2.2 (3) 




273.15 


17.62 


0.493 


(8) 


1.7 (2) 


0.178 (9) 


280 


17.46 


0.64 


(1) 


1.4 (2) 




298.15 


17.08 


1.07 


(1) 


1.1 (1) 


0.171 (9) 


308.15 


16.88 


1.33 


(2) 


1.0 (2) 




320 


16.66 


1.68 


(1) 


0.8 (1) 




328.15 


16.45 


2.01 


(2) 


0.61 (8) 


0.161 (8) 


333 


16.31 


2.20 


(2) 


0.60 (7) 




338.15 


16.26 


2.40 


(2) 


0.47 (7) 
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Table 4: Density, self-diffusion and Maxwell- Stefan diffusion coefficient as well 
as shear viscosity of the mixture methanol (1) + ethanol (2) at 298.15 K and 0.1 
MPa from present molecular dynamics simulation. The number in parenthesis 
indicates the statistical uncertainty in the last digit. 



XI p Di D2 Di2 ri 



al mol ^ 


mol 1-1 


10"9 


—2 —1 

m s 


lO-^m 2s 1 


lO-^m-^s-i 


10 3pa 





27.02 






1.07(1) 






0.143 


26.22 


1.37 


(2) 


1.12(1) 


1.4 (3) 


1.17(9) 


0.305 


25.63 


1.52 


(1) 


1.25(1) 


1.8 (3) 


0.76 (9) 


0.500 


25.14 


1.73 


(1) 


1.44(1) 


1.7 (2) 


0.62 (8) 


0.570 


24.81 


1.82 


(1) 


1.53(1) 


1.6 (3) 


0.64 (9) 


0.754 


24.51 


2.07 


(1) 


1.72 (2) 


2.3 (3) 


0.70 (9) 


0.796 


23.91 


2.14 


(2) 


1.82 (2) 


2.0 (3) 


0.59 (7) 


0.890 


23.58 


2.28 


(1) 


1.90 (2) 


2.6 (3) 


0.60 (8) 


1 


23.21 


4.95 


(3) 
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Table 5: Vapor-liquid equilibrium data for the mixture methanol (1) + ethanol 
(2) at 393.15 K and 433.15 K from present molecular dynamics simulation. The 
number in parenthesis indicates the statistical uncertainty in the last digit. 





P 




yi 




P' 




Ah^ 


mol mol^ 


MPa 


mol mol 


mol 1-1 


mol 1-1 


kJ mol-i 


r = 393.15 K 





0.551 (6) 







14.84 (2) 


0.200 (1) 


33.24 (7) 


0.150 


0.448 (9) 


0.207 


(5) 


15.62 (2) 


0.154 (3) 


33.86 (6) 


0.282 


0.482 (8) 


0.369 


(2) 


16.35 (2) 


0.168 (3) 


33.32 (7) 


0.492 


0.538 (9) 


0.589 


(2) 


17.54 (3) 


0.185 (3) 


32.67 (8) 


0.678 


0.588 (9) 


0.752 


(5) 


18.87 (3) 


0.212(3) 


31.92 (9) 


0.838 


0.572 (9) 


0.884 


(3) 


19.94 (3) 


0.205 (3) 


31.43 (8) 


1.0 


0.767 (5) 


1.0 




21.30 (3) 


0.294 (1) 


29.74 (8) 


r = 433.15 K 





1.32 


(3) 







13.58 (4) 


0.469 (2) 


28.41 (7) 


0.150 


1.29 


(2) 


0.192 


(4) 


14.34 (4) 


0.453 (7) 


28.52 (9) 


0.282 


1.38 


(2) 


0.339 


(5) 


14.96 (4) 


0.486 (7) 


28.08 (9) 


0.488 


1.48 


(1) 


0.549 


(5) 


16.09 (3) 


0.529 (5) 


27.50 (8) 


0.588 


1.54 


(2) 


0.659 


(5) 


18.87 (3) 


0.548 (7) 


27.28 (8) 


0.674 


1.58 


(2) 


0.734 


(5) 


17.22 (4) 


0.560 (7) 


27.07 (9) 


0.836 


1.64 


(2) 


0.872 


(2) 


18.31 (4) 


0.594 (7) 


26.55 (8) 


0.942 


1.72 


(2) 


0.957 


(1) 


19.15(5) 


0.633 (7) 


26.12(9) 


1.0 


1.89 


(3) 


1.0 




19.42 (4) 


0.711 (2) 


26.12(9) 
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Table 6: Volume, excess volume, enthalpy and excess enthalpy of the mixture 
methanol (1) + ethanol (2) at 298.15 K and 0.1 MPa from present molecular dy- 
namics simulation. The numbers in parenthesis indicate the statistical uncertainty 
in the last digits. 



xi V h 



ol mol ^ 


cm ^ 


mol ^ 


—3 

cm ^ 


mol ^ 


kJ mol 1 


J mol-i 





58.44 


(5) 






-43.46 


(4) 




0.070 


57.15 


(5) 


0.00 


(7) 


-43.20 


(4) 


3 (52) 


0.200 


54.94 


(5) 


0.01 


(6) 


-42.75 


(4) 


2 (49) 


0.270 


53.72 


(5) 


0.03 


(6) 


-42.50 


(4) 


-1 (47) 


0.350 


52.28 


(5) 


0.00 


(6) 


-42.22 


(4) 


-9 (46) 


0.420 


51.08 


(5) 


0.04 


(6) 


-41.95 


(4) 


20 (45) 


0.500 


49.65 


(4) 


0.02 


(5) 


-41.69 


(4) 


-8 (44) 


0.600 


47.91 


(4) 


0.04 


(5) 


-41.32 


(4) 


5 (44) 


0.697 


46.19 


(4) 


0.04 


(5) 


-40.96 


(4) 


20 (44) 


0.800 


44.34 


(4) 


0.00 


(5) 


-40.61 


(4) 


1 (45) 


0.905 


42.48 


(4) 


0.00 


(5) 


-40.24 


(4) 


-1 (47) 


1 


40.81 


(4) 






-39.90 


(3) 
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Figure 3: 
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